% The main use of this function is to Compute Matrix Q by wwb June 1st

function Q = findQ(rimS, p, Nmax)
p = [0.1,0.5,0.9];
p1 = p(1); p2 = p(2); p3 = p(3);
Q = zeros(8,Nmax+1);     % Should I add 1 here?
for i = 1:8
    Pre_i = sum(rimS(i,:));
    for m = 1 : Pre_i
        % Then, how to find the combination?
         for type1 = 0:rimS(i,1)
            for type2 = 0:rimS(i,2)
                for type3 = 0:rimS(i,3)
                    if type1+type2+type3 == m-1   % I use m-1 
                    	Q(i,m) = Q(i,m)+ ...
                            nchoosek(rimS(i,1),type1)*(p1^(type1))*((1-p1)^(rimS(i,1)-type1))*...
                            nchoosek(rimS(i,2),type2)*(p2^(type2))*((1-p2)^(rimS(i,2)-type2))*...
                            nchoosek(rimS(i,3),type3)*(p3^(type3))*((1-p3)^(rimS(i,3)-type3));
                    end
                end
            end
         end
    end
end
end